STOCHASTIC EXPONENTIAL INTEGRATORS FOR 
FINITE ELEMENT DISCRETIZATION OF SPDEs FOR 
MULTIPLICATIVE & ADDITIVE NOISE 



Gabriel J Lord*and Antoine Tambue ''' 
Department of Mathematics and Maxwell Institute, 

Heriot-Watt University, Edinburgh EH14 4AS, UK. 
January 14, 2013 



Abstract 

We consider the numerical approximation of a general second order semi-linear parabolic stochastic 
partial differential equation (SPDEs) driven by space-time noise, for multiplicative and additive noise. 
We examine convergence of exponential integrators for multiplicative and additive noise. We consider 
noise that is in trace class and give a convergence proof in the root mean square l? norm. We discretize in 
space with the finite element method and in our implementation we examine both the finite element and 
the finite volume methods. We present results for a linear reaction diffusion equation in two dimensions 
as well as a nonlinear example of two-dimensional stochastic advection diffusion reaction equation mo- 
tivated from realistic porous media flow. Parabolic stochastic partial differential equation, finite element, 
exponential integrators, strong numerical approximation, multiplicative noise, additive noise. 

1 Introduction 

We analyse the strong numerical approximation of an Ito stochastic partial differential equation defined in 
Q.<zW' . Boundary conditions on the domain £1 are typically Neumann, Dirichlet or some mixed conditions. 
We consider 

dX = {AX^F{X))dt+B{X)dW, X{Q)=Xo, t e[0,T], T>0 (1.1) 

in a Hilbert space H. Here A is the generator of an analytic semigroup S{t) :— e''^ ,t >0 not neccessary self 
adjoint. The functions F and B are nonlinear of X and the noise term W (x, f ) is a g- Wiener process defined 
on a filtered probability space (D,^,P, that is white in time. The noise can be represented as a 

series in the eigenfunctions of the covariance operator Q given by 

^ (1-2) 

where [qi,ei), i e N'' are the eigenvalues and eigenfunctions of the covariance operator Q and j3/ are in- 
dependent and identically distributed standard Brownian motions. Precise assumptions on A, F,B and W 
are given in Section|2]and under these type of technical assumptions it is well known, see IHIEIE] that the 
unique mild solution of ( |l.l| l is given by 

X{t)^S{t)Xo+ [' S{t~s)F{X{s))ds+ [' S{t^s)B{X{s))dW{s). (1.3) 
Jo Jo 
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Typical examples of the above type of equation are stochastic (advection) reaction diffusion equations 
arising, for example, in pattern formation in physics and mathematical biology. We illustrate our work 
with both a simple reaction diffusion equation where we can construct an exact solution 

dX = {V-I>VX~XX)dt + dW (1.4) 

as well as the stochastic advection reaction diffusion equation 

dX = (v DVX - V • {qX) - l^piY^ dt +XdW, (1.5) 

where D is the diffusion tensor, q is the Darcy velocity field |4) and A is a constant depending of the reaction 
function. The study of numerical solutions of SPDEs is an active research area and there is an extensive 
literature on numerical methods for SPDEs For temporal discretizations the linear implicit Euler 

scheme is often used ||251|29| . spatial discretizations are usually achieved with finite element |[30l [3111271 . 
finite difference |l25]|29l or spectral Galerkin method 1221 |5l |6l H [U . 

In the special case with additive noise, new schemes using linear functionals of the noise have recently 
been considered |l5]|6l|7]|8l[9l[T0l. The finite element method is used for the spatial discretization in ll9l [T0l 
and the spectral Galerkin in ||5j |6l |7] |8] . Our schemes here are based on using the finite element method 
(or finite volume method) for space discretization so that we gain the flexibility of theses methods to deal 
with complex boundary conditions and we can apply well developed techniques such as upwinding to deal 
with advection. One of our schemes is the non-diagonal version of the stochastic scheme presented in 
II22II32II and the other is the extension of the deterministic exponential time differencing of order one ll39l 
to stochastic exponential scheme. Comparing to the schemes presented in |9, 10| on additive noise, the 
results here are more general since the linear operator A does not need to be self adjoint and we do not 
need information about eigenvalues and eigenfunctions of the linear operator A. Furthermore we examine 
here convergence for Ito multiplicative noise for the exponential integrators, which has not so far been 
considered for SPDEs for these integrators. As in IJOJ , schemes presented here are based on exponential 
matrix computation, which is a notorious problem in numerical analysis ifTTl . However, new developments 
for both Leja points and Krylov subspace techniques fT2l [T3l [T4l [TSl [T6l [TTl have led to efficient methods 
for computing matrix exponentials. The convergence proof given below is similar to one in |28 | for a 
finite element discretization in space and backward Euler based method in time. The paper is organised as 
follows. In Section|2]we present the two numerical schemes based on the exponential integrators and our 
assumptions on jl.l) . We also present and comment on our convergence results. Section [3] contains the 
proofs of our convergence theorems. We conclude in Section|4]by presenting some simulations and discuss 
implementation of these methods. 



2 Mild solution, numerical schemes and main results 
2.1 The abstract setting and mild solution 

Let us start by presenting briefly the notation for the main function spaces and norms that we use in the 
paper We denote by || • j] the norm associated to the inner product (•, •) of the Hilbert space H. For a 
Banach space Y we denote by || .|| ^ the norm of 'f, L^f) the set of bounded hnear mapping from "Y loY 
and by L2(D, ^) the space defined by 

L2(]D), X^) := I V random variable with value in 1^: E||v||^^ = j \\v{(o)\\ydV{o)) . (2.1) 

Let 2 ://—>// be a trace class operator We introduce the spaces and notation we need to define the 
2" Wiener process. An operator T E L{H) is Hilbert-Schmidt if 

\\nm--= Eiir.,i|2<-, 
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where (ci) is an orthonormal basis in H. The sum in H-H^^ is independent of the choice of the orthonormal 
basis inH. Wedenotethe space of Hilbert-Schmidt operators from toi/byL^ ■■=HS{Q^I^{H),H) 

and the corresponding norm || .||^o by 

l|r 11,0 := WTQ'I^Whs = f E \\TQ"^eA , T eL^ 

ViGN'' / 

Let (jO be a process, we have the following equality using the Ito's isometry |[T| 

E|l f (pdWf^ f EMlods ^ f E\\(pQ'/YHsd^. te[0,T]. 
Jo Jo ^2 Jo 

Let us give some assumptions required both for the existence and uniqueness of the solution of equation 
( |L1| | and for our convergence proofs below. 

Assumption 2.1 The operator A is the generator of an analytic semigroup S{t) = e''^, f > 0. 

In the Banach space ^((— A)"/^), a e M, we use the notation ||(— A)"/^.|| =: ||.||cc We recall some basic 
properties of the semigroup S{t) generated by A. 

Proposition 2.2 [Smoothing properties of the semigroup [181] 

Let a > 0, j3 > and < 7 < 1, then there exist C > such that 

\\{~Afs{t)\\L(H) < Ct-P for t>0 

\\i-A)-y{i-smLiH) < ay for t>o. 

In addition, 

i-Afsit) = Sit)i~Af on S)i{-A)P) 

If P > Y then ^{{-A)^) C &{{-Ay), 
\\DlSit)v\\p < Cr'-(^-«)/2||v||„, f>0, vG^((-A)«/2) Z = 0,1, 

; d' 
where D, :— — r. 

dt' 

We describe now in detail the assumptions that we make on the nonlinear terms F,B and the noise W . 

Assumption 2.3 [Assumption on the drift term F] There exists a positive constant L > such that F is 
continuous in H and satisfies the following Lipschitz condition 

||F(Z)-F(}')|| <L||Z-}'|| V Z,YeH. 

As a consequence, there exists a constant C > such that 

||F(Z)|| < ||F(0)|| + ||F(Z)-F(0)||<||F(0)|1+L||Z||<C(1 + ||Z||) ZeH. 



Assumption 2.4 [Assumption on the noise and the diffusion term B] 

The covariance operator Q is in reace class i.e. Tr( Q) < °°, and there exists a positive constant L > such 
that B is continuous in H and satisfies the following condition 

||B(Z)-B(F)||^o <L||z-y|| yz,YeH. 
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As a consequence, there exists a constant C > such that 

||B(Z)||^o < ||B(0)||^o + ||B(Z)-B(0)||^o<||B(0)|l^o+L|lZ|l<C(l + |lZ|l) ZeH. 

Theorem 2.5 [Existence and uniqueness [D] 

Assume that the initial solution Xq is an pQ—measurable H~valued random variable and Assumption 
Assumption \2.4\ are satisfied. There exists a mild solution X to ( |l.l| l unique, up to equivalence among the 
processes, satisfying 

\\X{s)\fds<^y (2.2) 

For any p >2 there exists a constant C — C{p, T) > such that 

sup E||X(f)||''<C(l +£11X011"). (2.3) 

te[0.,T] 

For any p > 2 there exists a constant C\ — C\ {p, T) > such that 

E sup ||X(r)||''<Ci(l+E||Xor). (2.4) 

te[Q.T] 

The following theorem proves a regularity result of the mild solution X of 

Theore m 2.6 Assume that Assumption \2.3\ and Assumption \2.4\ hold. Let X be the mild solution of ( |7.7| 
given in (jli}. //Xq € L2(D, &{{-Afl^)), j3 G [0, 1) then for all t e [0, T],X{t) € L2(D, ^{{-A)P/^)) with 

1/2 , / ,\l/2 1/2^ 



(E||X(f)|!^) <c(i + (e||Xo||2) +(l+E||Xo|p 



^ 1/2 

Proof Recall that if X is the mild solution of ( 1.1 1, according to ( |2.1[ ) we need to estimate I E||X(f)|| 



and check that 

(E||X(Olir'^' 



Recall that the mild solution is given by 

X(f)=5(f)Xo+ [' Sit-s)F{X{s))ds+ [' S{t~s)B{X{s))dWis) 
Jo Jo 



then 

1/2 / „ _ „o\l/2 



(E||X(f)||2)'' < {E\\S{t)Xo\\jy\(E\\J^S{t-s)F{X{s))ds\\l 

= I + II + III. 
Since Xq e L2(D, ^((-A)'^/^)), j3 e [0, 1), we obviouly have 

/=(E||5(f)Xo||^)'^'<c(E||Xo||^)'^'. 



.1/2 




) + 


H 







As a consequence of Assumption |2 . 3 1 and the semigroup properties in Proposition 2.2 we have 

/ i-t \ 1/2 rt / \ i/l 

II=[E\\J^Sit-s)F{Xis))ds\\jj < J^^{E\\Sit^s)F{Xis))\\l^ ds 

< f (E\\i~A)P/'S{t^s)F{Xis))fy^'ds 



< C^'||(-A)'5/25(f-.)||,„,(^„^.) (^sup^(Eil + \\Xis)\\ff'^ 

< c( ['it-s)-^ds] (l + sup (E||X(i)||2)'/^ 



< C 1 + sup {E\\X{s)\ 

0<s<t 



Q<s<t 
2\l/2' 
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Finally, Ito's isometry and the consequence of Assumption |2.4| yields 

II f = E|| [' S{t-s)B{X{s))dW{s)\\l 



E\\{-Af/^Sit~s)B{X{sm^.ods 

^2 

< c(^l^\\{~A)P/^S{t-s)\\l^^,^^^^ds^ (^1+ mpE||X(.)||2^ 



thus 



< C( 1+ sup E\\X{s) 

0<s<t 



l/2^ 

///<C I 1+ ( sup E\\X{s)"^ 



\Q<s<t 

Then, using ( |2.3| ) with p — 2 yields 

E\\Xit)\\jf' < c(i + (e||Xo|12)'^'+ mp^(E||X(.)||2)'/' 



< C (^l + (^E||Xo||i) +(l+E||Xo| 

<c 

Then, if Xq G L2(D, ^((-A)^/^)), j3 € [0,1) then for all f € [0,T],X{t) € L2iI]>,9{{-A)Pf^)). 
More results about the regularity of the mild solution X can be found in BOl 1411 . 

2.2 Application to the second order semi-linear parabolic SPDEs 

We assume that Q. has a smooth boundary or is a convex polygon of W' , d = 1,2,3. In the sequel of 
this paper, for convenience of presentation, we take A to be a second order operator as this simplifies the 
convergence proof. 

More precisely we take H (H) and consider the general second order semi-linear parabolic stochas- 
tic partial differential equation given by 

dX{t,x) = {V-DVX{t,x)-q-VX{t,x)+f{x,X{t,x)))dt + b{x,X{t,x))dW{t,x), (2.5) 

x£Q.,t E [Q,T] where f,b:Q.xR^Raie two continuously differentiable functions with globally bounded 
derivatives. 

In the abstract form given in the linear operator is defined by 

where we assume that Di j €E L°° (H) , qi € L°° (i2) and that there exists a positive constant ci > such that 

d _ 

£ Dij{x)^i^f > ci 1^ |2 V<^ e xeQ. ci > 0, (2.7) 

'■.i=i 

andF:H^H,B:H^ HS{Q^I^{H),H) are defined by 

{F{v)){x)^f{x,v{x)), {B{v)u){x)^b{x,v{x))-u{x), (2.8) 

for all jc e £2, v e //, m G Q^^^{II), with H = L^(£2). As the functions / and b are two continuously differ- 
entiable functions with globally bounded derivatives, the Nemytskii operator F corresponding to / and the 
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multiplication operator B defined in (|2.8|l satisfy Assumptions |2 . 3fp4| for appropriate eigenfunctions such 



that sup 



sup|e,-(x)| 



< oo (see |40, Section 4 ]). 



Notice that by the definitions of the operator B and 1 1 . 1 1 , for Y £ H = L^{Q.) 

\\B{Y)\\l,= '£\\b{Y)Q'/^e.f, (2.9) 

where b{Y) is the Nemytskii operator defined by 

b{Y){x)^b{x,Y{x)) xGH. (2.10) 

We introduce two spaces H and V where H C V that depend on the choice of the boundary conditions 
for the domain of the operator A and the corresponding biUnear form. For Dirichlet boundary conditions 
we let 

V^M^H^ (Q.) = {v e //' : V = on dQ.}, 
and for Robin boundary conditions, Neumann boundary being a special case, we takeV =//'(i2) and 

M={veH\Q.):dv/dvA + aov = on dQ.} , Oq G M. 
See 1191 for details. The corresponding bilinear form of —A is given by 

f f 4- du dv 4, du \ 

a[u,v) = Vl^qi-^vXdx u,v eV (2.11) 

Jn oxj oxi dxj J 

for Dirichlet and Neumann boundary conditions, and by 

, , f ( 4^ du dv 4^ du \ f 

a{u,v) ^ / ^Dij- — ^ y^cfi^ — v\dx+ / oCQUvdx u,v£V, (2.12) 

Jn ' oxj dxi dxj I Jdn 

for Robin boundary conditions. According to Garding's inequality (see ll23l[T9l ). there exists two positive 
constants cq and Aq such that 

a{v,v)+co\\vf>XoMliia) ^ ^- (2.13) 

By adding and subtracting coXdt on the right hand side of ( |l.l| i, we have a new operator that we still 
call A corresponding to the new bilinear form that we still call a such that the following coercivity property 
holds 

fl(v,v)> Ao|lv||^,(^) Vvey. (2.14) 



Note that the expression of the nonlinear term F has changed as we include the term —cqX in a new 
nonlinear term that we still denote by F . The coercivity property ( 2.14| i implies that A is a sectorial on 



L?{Q.) i.e. there exists Ci, 9 £ {hn,n) such that 



||(A/-A)->|l^(^2(n))<^ (2.15) 

where 5e = {A e C : A = pe"t', p > 0, < |^| < 0} (see |[H[19]). Then A is the infinitesimal generator 
of bounded analytic semigroups S{t) : = e''^ on L^(£2) such that 

S(t):=e"^ = J- [ e'^(Xl-A)-^dX, f > 0, (2.16) 
2 71;; j'^f 

where denotes a path that surrounds the spectrum of A. 

Functions in H satisfy the boundary conditions and with H in hand we can characterize the domain of 
the operator (— A)*^/^ and have the following norm equivalence If36l[37ll for r = 1,2 

||y||^,-(a) ^ ||(-A)'-/2v|| ^: ||v|U Vv G f^((-A)'-/2) = Hn//'-(il). 
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2.3 Numerical schemes 

We consider the discretization of the spatial domain by a finite element triangulation. Let ii^, be a set of 
disjoint intervals of Q. (for d = I), a triangulation of Q. (for c/ = 2) or a set of tetrahedra (for d — 3) with 
maximal length h. Let Vi, CV denote the space of continuous functions that are piecewise linear over the 
triangulation J^,. To discretize in space we introduce the projection f/, from L^(i2) onto V), defined for 

u£L^{D.) by 

{Phii,x) = iu,X) yX&Vh. (2.17) 
The discrete operator A/, : V/, V/, is defined by 

{Ah(p,X)^{A(p,X)^-a((p,X) (P:X(^Vh. (2.18) 

Like the operator A, the discrete operator A/, is also th e gen erator of an analytic semigroup Si, ■= e'^'' ■ 

The semi-discrete in space version of the problem ( 1.1 1 is to find the process =X^{.,t) e V/j such 
that for f e [0,r], 

dX*" = {AhX^ +PhF{X^))dt +PhB{x'')dW, X''{0) = PhXo. (2.19) 
The mild solution of ( 2.19[ ) at time t,„ = inAt, Af > is given by 

X''it,„) = Sh{t,n)PhXQ+ f'" Sh{t,n-s)PhF{X\s))ds 

Jq 

Sh{tm-s)PhB{x'')dW{s). (2.20) 





Then, given the mild solution at the time t,„, we can construct the corresponding solution at f„,+i as 

X''{t,n+i) = Sh{At)x'\t,„)+ Sk{At-s)PhF{x''is + t„,))ds 

Jo 

Sh{t,n+l~s)PhB[X'')dW{s). 

II 

To build the first numerical scheme, we use the following approximations 

ShiAt-s)Fix\t,„+s)) « Sh{At)F{x''{t,„)) se[0,At], 
Si,{t,„+i-s)P^B{x''{s)) « Sh{At)P,,B{x''{t,n)) se[t„,t„+i]. 

We can define our approximation Y/^ ofX{mAt) by 



= (po(AfA,,) [Y;;,+PhF{Y^,) +P,B(y„';)AW,„ ) , (2.21) 

where 

(poiAtAi,) := 



AW,n := W,„+i - W,n = VAf V^iRi,mei, 

with Ri m are independent, standard normally distributed random variables with means and variance 1 . 
We call the scheme in ( |2.21[ ) SETDMO. To build the second numerical scheme, we use the following 
approximations 

F{x\t„+s)) « F{x'\t„)) .€[0, Af], 
Si,itm+i-s)PhB{x''{s)) w Si,{At)PhB{x\t„,)) s e [t,n,t,n+l]. 
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We can define our approximation ofX{mAt) by 

= /'^"Xi; +A/;' (/'^" -/) PhF{xi)+e^''"PhB{^i) - . (2.22) 

For efficiency we can rewrite the scheme ( |2.22 1 as 

^m+i = x!'„+At(p, (AtAh) (Ah (xl;,+PhBix!:,)Aw,„) +PhF{x';„)\ 

where 

(piiAtAh) = {AtAh)-^ (e^^K -/) = 

We call the scheme in ( |2.22[ ) SETDMl. This scheme is also used in |24| with the Fourier method to solve 
fourth order stochastic problems. 

2.4 Main result 

Throughout the paper we take f,„ = mAt S (0, T], where T = MAt for m,M G N. We take C to be a constant 
that may depend on T and other parameters but not on At or h. Our result is a strong convergence result in 
l2 for schemes SETDMl and SETDMO. 

Theo rem 2.7 Let X(t,„) be the mild solution of equation \\..\) at time t,„ ~ mAt, Af > represented by 



{1.31 Let l^l^ be the numerical approximations through ( |2.22| l or \2.2\) (Cm— scheme SETDMl and 



Q^^YllJor scheme SETDMO) andQ<y<\. Assume f/zaf Z7(L2(D, ^((-A)"))) C L2(D, ^((-A)"))/or 
(X £ (0, 7/10) small enough. The following estimates hold. 
//XoeL2(D,^((-A)1')), 0<7< 1/2 then 

{E\\xit,„)-Cfy' < c(fi-'+^^)/^/.+Afr/2). 

//XoeL2(D,^((-A)1')), l/2<7<lf/zen 

(E\\x{t,„)-Cfy^' < c[h+Aty/^). 

Suppose f/iafF(L2(D,^((-A)"'))) CL2(D,^((-A)"')), a' E (0,7/10) small enough: 
IfXo e L2(D,^((-A)^)) then 

{E\\x{t,„)-Cff' < c{ti-'^'^h' + Atr/'). 

IfXoeL2(D,&{-A)) andb{L2{D,S'{-A))) CL2{D,^{-A)) then 

1/2 



E\\X{tm)-Cf) < c(/i2+Af(i/2-^)). 



£ G (0, 1 /2), small enough. 

Remark 2.8 In the proof of Theorem^ the assumptions b{L2{D, S^{{-A)"))) C L2(D,^{{-A)'^)) for 
a e (0,7/10) andF{L2{D,^{{-A)"'))) C ^2(0, ^((-A)"')), a' & (0,7/10) are enough. We only need 
that b(X{t)) e L20Di,&{{-A)")), a e (0,7/10) andF{X{t)) e L2(D, ^((-A)"')), a' G (0,7/10), Vf G 
[0, T], where X is the mild solution of equation ( |l.l| l. 

Although we have taken the linear operator A to be a second order operator, similar results will hold, 
for higher order operators. Computationaly, the noise given by \1.2) is truncated to terms. Therefore 
the corresponding approximated solutions become Xm'^ for SETDMl and Y^ '^ for scheme SETDMO. For 
noise where the eigenvalues of the covariance operator have a strong exponential decay, Xm'^ and Fm'^ 
are close to X^' and Fjj respective ly. I n the case of additive noise, it has been proved in f3TJ that with the 
truncation to terms of the noise d 1 .2 1 the corresponding discrete mild solution X''-^ in ( 2.20 1 has the same 



order of accuracy respect to h as X^ 

We also note that smooth noise improves the accuracy in Theorem 2.7 for additive noise (see |32 | and 
Figure [Tjin Section|4]|. 
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3 Proofs of the main results 
3.1 Some preparatory results 

We introduce the Riesz representation operator Rh :V ^ Vj, defined by 

{-ARhV,x) = i-Av,x) = a{v,x) v £ V, V;i: G V,,. 

Under the regularity assumptions on the triangulation and in view of V— elHpticity ( |2.7[ ), it is well known 
(see lfT9l [38 n that the following error bounds holds 

||/;/,v-v||+/;||/;,,v-v||^i(jj<C/;n|v||;,,.(n), veVnH'-{£l), re {1,2}. (3.1) 

We start by examining the deterministic linear problem. Find u eV such that such that 

u'^Au given m(0) = v, t e {Q,T]. (3.2) 

The corresponding semi-discretization in space is : Find m/, e V/, such that 

M/, — AhUh 

where = f/,y. Define the operator 

r,,(f) := S{t) - Sh{t)Ph = e'^ - e'^"Ph (3.3) 

so that u{t) — Uf,{t) ~ Th{t)v. 

Lemma 3.1 The following estimate holds on the semi-discrete approximation o/( |3.2[ ). Ifv e S>{{~A)^^^) 

\\u{t)-uh{t)\\ = ||r,(r)v||<C/2'-f-('-/^)/2||y||^ re{l,2}/3<r, (3.4) 

where r is related to ( |3.1| l. 
Proof 

The proof for r = 2 and P — can be found in A^TPI/. Theorem 7.1, page 817]. 

Set 

Uh(t)~u(t) = {uh{t) ~ Rhu{t)) + {Rhu{t) - u{t)) = e{t)+p{t). (3.5) 

It is well known [36 J that A^Rh ~ PhA. Indeed for v G S!{A), % G have 

{PhAv,x) = {Av,x) (by definition of Ph) 

= {ARhV.X) {by definition of Rh) 

= {AhRhV, X) [since RhV e V/,) 

thus A^Rh = PhA. We therefore have the following equation in 9 

e,=A,,d-PhD,p. 

Hence 

d{t)^Sh{t)d{Q)- f Sh{t-s)PhDspds. 
Jo 

Splitting the integral up into two intervals and integration by parts over the first interval yields 
B{t)=Sh{t)d{Q)+Sh{t)Php{0)-Sh{t/2)Php{t/2)+ f'^ {DMt-s))Php{s)ds- f Sh{t ^ s)PhD,p{s)ds, 

Jo Jt/2 
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with Ds = d / ds. Since 9{t) (EVh we therefore have Phd{t) ~ 9 (f ), then 

rt/2 



B{t)^Sh{t)PhTh{0)v-Sh{t/2)Php{t/2)+ f {DMt~s))Php{s)ds- f Sh{t - s)PhD,p{s)ds. 

Jo Jr/2 



Since 

nr,,(o)v-p,,(v-p,,v) = o, 

we therefore have 

d{t) = -Sh{t/2)Php{t/2)+ f'^ D,Sh{t-s)Php{s)ds- f Sh{t-s)PhD,p{s)ds. 

Jo Jt/2 

Using the fact that Si, and Pf, are uniformly bounded independently of h with the smoothing property of Sh 
in Proposition \2.2\ yields 

110(011 <c(^||p(r/2)||+^'^'(f-.)-ipWI|fl'. + ^^J|D,pWII«'^' 



The estimate with the smoothing property ofS{t) in Proposition 2.2 yields 
||p(f)|| < C/i1m||, < C/i'"f-('-^)/2||v||^ 

\\D,p{t)\\<Ch'^\DM\r<Ch'-t-'-^'--S')l^v\\p, re {1,2}, /3 < r, v G ^((-A)'5/2). 



Then 



\\e{t)\\<ch'-t-('--py^v\\p+ch'-\\v\\py^^\t~s^^^ 

{t-s)-'s-^'-py'ds+ f s-'-^'-py^ds<ct-('-py\ 

Jt/2 



Since 

rt/2 



t/2 

we therefore have 

l|7},(Ov'll<l|0(OII + IIP(OII<c/.'-r('-'5)/2||,||^ 

Our second preliminary lemma concerns the mild solution SPDE of ( |1.1[ ). 



Lemma 3.2 Let X be the mild solution given in {1.31 Suppose that Assumption 2.3 and Assumption 2.4 



hold. Let < Y < I, fi,f2 G [0,r] be so that ti < t2. IfXo e L2{V),&{{-Ay)) then we have the following 
estimate, 

nnt2)^x{t,)f < c{t2-t,y (e\\xo\\i+e( sup (i+iix(.)ii)y+( sup E(i+iix(.)iin ). 

Y V0<i<7' / \0<s<ti I I 

Proof Consider the difference 

X[t2)-X{t{) 

= (5(f2)-5(fi))Xo+ S{t2-s)FiX{s))ds- j\iti-s)F{X{s))ds 

+ (^j'^ S{t2 - s)B{X)dW{s) - j^'' S{ti - s)B{X)dW{s) 
= I + II + III, 

so that E\\X{t2) -X{ti)\\^ < 3(E||/||2+E||//|p + E||///||2). We estimate each of the terms /,// and ///. 
Estimation of the terms / and // are similar to ones in [||9i ^lOJ, Lemma 3.2] with additive noise. Using 
Proposition 2.2 as in ||9l[T0l yields 



||/|| = ||5(fi)(-A)-r/2(l-5(f2-fi))(-A)r/2xo|| < C{t2^t,y/^\\Xo\\ 



7- 
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and 

E||//|l2<C(f2-fi)'^E( sup (1 + ||X(^)||)" 



Vo<.s<r 



For term ///, we have 

iSit2-s)-Siti-s))B{X)dWis)+ r 

10 Jti 

Using the Ito isometry property yields 



/// = /'' {S{t2 - s) - S{ti - s))B{X)dW{s) + r S{t2 - s)B{X)dW{s) = IIIi +III2. 

Jo Ju 



E||///i||2 = E\\ r {Sit2-s)-S{ti^s))BiX)dW{s)\\ 







'I 



^EUS{t2-s)-Sin-s))BiX)Q'/YHsd^- 
Using Assumption 2.4 and Proposition |2 . 2| yields 







E||///i||2 < Cy^'\\{S{t2-s)-S{t,-sml^,^^^^ds^ ^up E(l + ||X(.)||)2j 

= C^"||5(fi-.)(-A)''/2(-A)-^/2(I-5(f2-fi))|l^(^2(^))«'.) (^^sup E(1 + ||X(.)||)2 
= C^"||(-A)r/25(fi-.)(-A)-r/2(I-5(f2-fi))||2(^,(^))J.) ( sup E{l + \\X{s)\\f 



, 0<.s<f| 



< c{t2-hy(^J\n~s)-^ds^ 



sup E(1 + ||X(.)||)^ 

.Q<s<ti 



< C{t2-hy{ sup Eil + \\Xis)\\Y 

\ 0<s<?i 



Let us estimate E||///2||. The Ito isometry again, with the boundedness of S and Assumption |2.4| yields 



E||///2||2 = E\\ ["s{t2-smx)dW{s)f 

"E\\Sit2-s))BiX{sml2ds 
H 



r||5(f2-5)||^i2(n))«'^( sup E{l + \\X{s)\\f 

2 



< C{t2-h) \ sup E(l + ||X(i)||) 

yo<.v<fi 



Hence 



E||///|l2<2(E|l///i|l2+E|l///2|l2)<C(?2-fi)^. 
Combining our estimates of E||/|p,E||//|p and E||///|p ends the lemma. 

3.2 Proof of Theorem |2J] for the scheme SETDMl 
Proof Set 

rtk+i ft,,, 
X{t„,) = S{t,„)Xo+Y, Sit,„-s)F{X{s))ds+ S{t,„~s)B{X{t,n))dW{s) 
k=o ■''k Jo 

= X{t^)+0{t,n). 
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Recall that 

= r e^^'-"^^>-PhF{Xi_,)ds+ ['"' e^''"-'^'^^PhB{X,'l_i)dW{s) 

Jo Jt„,_i 

= Sh{t,n)PhXo+Y, { Sh{t,n-s)PhF{xl')ds+ Sh{t,„-s)P,,B{xl')dW{s) 

"1-1 / ! \ "1-1 rti^^i 

= Sh{tm)PhXa+Y^( Sh{tm~s)Pt,F{4)ds)+Y. Sh{t,n'S)PhB{Xt)dW{s) 

= 7'' + o'^ 



with 



"1-1 / 1 

= Sh{t,„)PhXQ+Y,{ Sh{t„~s)PhF{xl')ds 

k=Q y-'tk 



We examine the error 



x{t,„)^x';„ = xit„,)+oit,„)-x;;, 

= x{t„,) + o{t,„)-(zi:^+ot) 

= (x{t^)-z':,) + [oit,„)-o',i 

= I + 11, (3.6) 



thus 



nntm) -Kf < 2 (Eii/||2 +Eii//|p) . (3.7) 

We follow the approach in [10|. Let us estimate the first term E||/|p. Using the definition of T), from 
|3.3|l, the first term / can be expanded 



"1-1 rtk^[ 

I = ThXo+Y, S{t„,-s)F{X{s))-Sh{t,„-s)PhF{x'[)ds 



k=0 -"k' 
"-1 ftk+l 



"'—1 nk+l 

ThXo+Y, / Sh{t,„-s)Ph{F{X{tk))~F{X^^)))ds 
+ E / Sh(t„,~s)Ph{F{X{s))-F{X{tk)))ds 

k=0 •''k 



Then 



+ ^ / {S{t,„^s)~S,,{t,„-s)P,,)F{Xis))ds 

k=Q ■''k 

= h+h+h+h- (3.8) 



E||/||2<4(E||/i||2 + E||/2||2+E||/3||2+E||/4|| 



Let us estimate /i, for < 7< 1 with 27 < r and r e {1,2}, if G L2(D, ^((-A)^')), Lemma [iT] with 
j3 = 27 yields 
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and if Xq e ^{—A) we have 



Epilp <C/!%||Xop 



2- 



For I2, using Assumption 2.3 triangle inequality as well as the fact that Sh{t) and P/, are bounded 
operators with Fubini's theorem yields 



EII/2II- < Cm J2 E|| / 5/,(f,„-^)PjF(X(f,))-F(X,'') W^|| 

k=0 •''k ^ ' 

<r=0 V^'t 



< CmAf £ / ^ (E||X(r,)-Xi'||2)flf, 



m— 1 

Once again using the Lipschitz condition, triangle inequality, the fact that 5/, and P/, are bounded but 
with Lemma [J!2] yields 



m-l 



m— 1 



< 



k=Q -"k 



\k=o-''k J 



thus 



(2 \ '/2 

E||XoE + Ef sup (l + ||X(i)||)) +f sup E(1 + ||XM||)2' 
Vo<.s<7' / \o<.v<r 

< CAt^/^(E\\Xo\\l + E( sup il + \\X{s)\\)) +f sup E(1 + ||X(.)||)2' 

y Vo<i<7' / Vo<i<7' 



E\\hf<CAt^. 

r— 1 ^ ^ 



If € L2 (D,^(-A)) we obviously have Epalp < C(Af)'"'^ by taking 7= 1 - e in Lemma 3.2 
(0, 1 /2) small enough. Let us estimate (E||/4|p) . For r = 1,)3 = 0, using Lemma 3.1 yields 



thus 



(E||/4f)'/' < f""\E\m^-s)FiX{sm'f'ds 

k=Q •''k 

< Ch sup {E\\F{X{s))ff'( f'"'{t,,~s)-'/^ 

0<s<T \J0 

< Ch, 



E\\hf <Ch^. 
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If F(L2(D,^((-A)« ))) c L2(D,^((-A)" )), a' E (0,7/10) small enough, for r = 2,j3 = 2a', using 
Lemma |3.1| yields 



1 /? f'k+l / \ 1/2 

(E||/4||2)'/' < Y. / (E||7),(f„,-.)^'(X(.))||2) ds 

< Ch^^sup^(E\\F{X{s))\\}y'^ (^fj {t„,-s)-'+l^l^ds^ 



< Ch^, 

thus 

2 ^ r-ji 



E||/4ir <c/!" 

Combining the previous estimates yields: For Xq E L2 (D, f^((— A)'')) , 1/2 < 7 < 1 



E 



ForXoeL2(D,^((-A)'')), 0<7< 1/2 

For e L2 (D, Si{{-Ay)) , < 7 < 1 and if F(L2(D, ^((-A)«'))) c ^2(0, ^((-A)«')), a' e (0, 7/IO) 
small enough 

/ m-l ,.t,. , , , 

^/;||2 



f„2+2r;,4^^^r+^ / U\\x{t,)-xt 

ForXo eL2(D,^(-A)) and if f (L2(D, ^((-A)«'))) c L2(D, ^((-A)"')), a' E (0,7/10) small enough, 

/;4 + Af'-^+^ {E\\X{t,)-Xj:f)ds\, 

with e e (0, 1 /2) small enough. 

Let us estimate E||//|p, we follow the same approach as in 1281 . Note that in the case of additive 
noise the estimation is straightforward and smooth noise improve the accuracy (see f32l and Figure [T] in 
Section|4|. For multiplicative noise we have 



r'm "1-1 ftk+i , 

// = / S{t,„-s)B{Xis))dW{s)-Y, Sh{t,n-tk)PhB{xi)dW{s) 

Jo Jrt 

= L / Sh{t^-h)pJB{X{tk))-B{xj;)))dW{s) 
+ E / Sh{t,,-h)Ph{B(X[s))-B{X{tk)))dW{s) 

k=0 •''k 

+ ^ / {S{t„-ti,)-Sh{t,n-tk)Ph)B{X{s))dW{s) 

k=0 •''k 

+ ^ / {S{t„,^s)^S{t^^tk))B(X{s)))dW{s) 

k=0 ■''k 

= Hi +112+113 +114. (3.9) 
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Then 

E|l//|l2<4(E||//i|l2+E||//2|l2+E|l//3||2+E||//4|l2). 
Let us estimate each term. Using the Ito isometry, the boundedness of 5/, and P/,, and Assumption |2.4| yields 

r'k+i 



E||//i||2 = E\\J^ r Shit„,~t,)P,,(BiX{t|,))^Bixj:))dW{s)\\■ 
^E|| y 5/,(f„-r,)i';,(B(X(f,))-B(xi'))^W(^)||- 



< C£ / E||B(X(f,))-B(Z,^)|L2«f. 



m— 1 



< C£ / E|lX(f,)-X/'||2flf,. 

A:=0 



For € L2 (D, ^ ((-A)^)), using Lemma 3.2 and Assumption|Z4|yields 



n-l 



E||//2||2 = E||^ / Sh{U-h)Ph{B(^X{s))^B{X{tk)))(m{s)\\ 

k=0 



= £ / E||5/,(f„,-r,)n(B(X(^))-B(X(f,)))||^ofl'* 

"'-1 ! 

< / E||X(^)-X(fi)!|2«fs 

A:=0 



' ni — 1 



< c £ / (5-f,)^rf. 

X (e\\Xo\\1 + E ( sup (1 + ||X(^)||)) + ( sup E|| (1 +Xis)\\y 

\ \0<s<T J \q<s<T 

< CAtTf. 

For Xq g L2 (D, ^(— a)), taking 7 = 1 — e, with e small enough yields 

E||//2||^<CAf'"''. 
Let us estimate E||//3|p. By Ito's isometry and Lemma 3.1 we have 

n'-l filial 

E\\Ihf - nL {S{t„-tk)~Sh{t,„-tk)Ph)B{X{s))dW{s)f 

= / ms{tm-tk)-Sh{t^-tk)Ph)B{xm\i2ds 



1 



£ / E||7},(f,„-fi)B(X(i))||22£!!i. 



Indeed using Lemma[3l]andAssumption[24| for fc(L2(D,^((-A)")))cL2(D,^((-A)")), ae (0,7/10) 
small enough, we have 

\mt,n~h)Bixismi2 = f;iir,(f,„-f,)fe(x(i))ei/2e,-i|2 

" i€N 

< Y.\\T^Hit,n-mX{s))f\\Q'/^e£ 

< Ch^{t„,-t,)-'+"\\b{Xis))\\lTriQ), 
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thus 

/ m-l ^ 

E\\Ih\\^ <Ch^Tr{Q) sup E\\b{X{s))\\l Af« ^ (m - /t)-i+" 

0<s<T \ k=o J 

since 

m— 1 

is the discrete form of 

f-m— 1 

Jo 

we therefore have 



niW<ChHAQ) sup E|lK^(5))|l2. 

0<s<7' 



For Z7(L2(D,^((-A)))) C L2(D,^((-A))), we obviously have using Lemma 3.1 



E||//3||2<C/i4Tr(e) sup E||K^(i))||i 



0<.v<7' 

2 



Let us estimate E||//4|| , by our assumption on h the following estimation holds 



m— 1 



E||//4f = E / n\{^{Un~s)~S{t„,~h))B{X{s))f^^ds 

k=0 ■''k " 



E / \\{-A)-'^l\S{t,,,-s)-S{t^-ml , d^^ sup E|1KXW)|1 



2 

a I I 



smce 



\\i^A)-"/^{S{t„,-s)-Sit„,-t,m%L2m = ll(-A)(>-«)/^5(f„,-.)(-A)(-'/2)(I-5(,v-f,))ll^(i2(n)) 

< Cis~tk)it„,-sY"-'\ 



thus 



E / (5-f,)(f„~^)(«-')rf. sup E\\b{X{s))\\l 

k=0-'tk J \0<s<T 



< At[ supE\\biX{s))ra\). 

\0<s<T / 

Combining the estimates related to // yields the following. 

That for Xq G L2 (D,D(-A)1') and b{L2iD, Si{{~A)"))) c L2(D, ^((-A)«)), a e (0,7/10) small 
enough, 



n-l 



tk+l 



E\\II\\^<C[h^+Atr+Y, / E\\X{t,)^X^\\^ds]. 

\ k=0 •''k J 

ForXo €L2(D,D(-A)'') and /7(L2 (D, ^(-A))) C L2 (D, ^(-A)) 

h^+Aty+ E / E\\X{tk)-X'i\\^ds\. 

k=0 •''k J 

Combining the estimates of E||/|p and E||//|p and applying the discrete Gronwall lemma ends the proof. 
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3.3 Proof of Theorem |2J] for the scheme SETDMO 

We just give a sketch of the mains steps. Recall that 

JO J Tff,— \ 

= 5/,(f,„)P/,Xo+£ / Sh{t,„-tk)PhF{Yl')ds+ Sh{tn.-tk)PhB{Yj:)dW{s) 

'"-1 / /-rj-^i \ m-1 /-fj.+ i 

= 5,,(f„,)P/,Xo+E / Shit,„-tk)PhF{Y^^)ds) + Y, Sh{t,n-tk)PhB{Yl')dW{s) 
= z''+o''. 



We can therefore put the estimation of the error in the form of equation |3.6| The estimate of the corre- 
sponding E||/|p is the same as in Theorem : 



2.7 



with the extra term 



tk+l 



m— 1 

^5 - E / " " (^(f™ -')- ^(^'« - tk))F{X{s))ds. 

k=0 •''k 

This is estimated in ||9l Theorem 2.6 as 

(Epsll^)'''^ <C(Af + Af|ln(Af)|) <CAf1'/2. 
The estimation of E||//||^ is the same as for the scheme SETDMO. 

4 Simulations 

Efficient implementation of (p,, / = 0, 1 can be achieved by either the real fast Leja points technique in ifTSl 
[T6l [TTl [TOl or the Krylov subspace technique in lfT2l [13] [TOl . In the first example we apply the scheme 
to linear problem where we can construct the exact solution for the truncated noise. The finite element 
method is used for space discretization. In this example we use the real fast Leja point technique to 
compute the exponential functions (pi, / = 0, 1. We use noise with exponential correlation (see below) 
which is obviously a trace class noise. In the second example we apply the scheme to nonlinear stochastic 
flow with multiplicative noise in a heterogeneous media. To deal with high Peclet number flow, we use the 
finite volume method for the space discretization. In this case we use the Krylov subspace technique to 
compute the exponential functions (pi, / 0, 1, implemented in the matlab functions expv . m and phiv . m 
of the package Expokit [13|. We compute the exponential matrix functions 9, with the Krylov subspace 
technique with dimension m ~ 6 and the absolute tolerance 10^^. 

In the legends of our graphs, "SETDMl" denotes results from the SETDMl scheme, "SETDMO " 
denotes results from the SETDMO scheme with "Implicit" denotes results from the standard semi-implicit 
Euler-Maruyama scheme. 

As a simple example consider the reaction diffusion equation with additive noise in the time interval 
[0,T] with diffusion coefficient D > 

dX = {DAX-Q.5X)dt + dW X{0)^Xo, £2 = [0,Li] x [0,L2] 

with homogeneous Neumann boundary condition. As the exact solution is known for comparison, we take 
/ in the equation (|2.5[) to be linear here 



/(m) = -0.5m. (4.1) 
The corresponding Nemytskii operator F is obtained from p.8|l. Of course, in general, F will be nonlinear. 



F verifies obviously Assumption 2.3 Here b{x,u) — I, x € Q., u € 
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We consider the covariance operator Q with the following covariance function (kernel) which has strong 
exponential decay 



Cr{{xi,yi);ix2,y2)) 



Abib2 



exp 



{x2-xi) {y2~yi) 



where b\,b2 are spatial correlation lengths in x— axis and y- axis respectively and F > 0. 
It is well known that the eigenfunctions {e|'^e^^^'},j>o of the operator A = DA is given by 



^0 



J (0 
A,) 



0, e 



(') 



2 iiC) \ 
— cos(A> x) 



I7t 

Li 



/e{l,2} /=1,2,3,--- 
with the corresponding eigenvalues {Xij}ij>o given by 



The corresponding values of {li,j}i^j^Q in the representation (1.2 1 are given by 

Fexp 



L((A(').OH(AfV) 



(4.2) 



2n 

see Il23l for details and 1251 l26l . We compute the exponential functions (p,, / = 0, 1 with the real fast Leja 
point technique and the absolute tolerance 10^^. In our simulation we take Li — L2 — I and the finite 
element triangulation is contructed with the rectangular grid with size Ax — Ay=l/ 150. Figure 1(a) shows 
the time convergence of SETDM1,SETDM0 and semi-implict schemes. The three methods have the same 
order of accuracy. The temporal order of convergence that we observe is 0.9 for all the schemes. This is 
higher than the predicted theoretical order of comvergence 0.5 in Theorem |2.7| The noise is regular and 
this order agrees with that in \ 




log (At) 



(a) (b) 

Figure 1 : (a) Convergence of the root mean square L^ norm at T = 1 as a function of At with 10 realizations 
with Xo = 0,F=l,D=l. The noise is white in time and with exponential correlation in space with lengths 
b\ =b2= 0.2. The temporal order of convergence in time is 0.9 for all schemes. In (b) we plot a sample 
true solution. 

As a more challenging example we consider the stochastic advection diffusion reaction SPDE with 
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multiplicative noise 

dX = ( V-DVX- V- (qX)- 



X 



X +1 



D 



10-2 






10-3 



dt+XdW, Xo = a = [0,1] X [0,1] (4.3) 



(4.4) 



with mixed Neumann-Dirichlet boundary conditions. The Dirichlet boundary condition is X = 1 at .jc = 
and we use the homogeneous Neumann boundary conditions elsewhere. According to Theorem 2.5 we 
need to take the initial data 

XoeL2(D,^{{~A)l^)), j3 >Otohave a regular solution such that XdW make sense. For our simulation 
we take Xq = 0. For a homogeneous medium, we use the constant velocity q = (1,0). In terms of equation 
(|2.5|l the nonlinear terms / and b are given by 



f{x,u) 



u 

\u\ + l) 



b{x,u) — u, M G M, X G £2, 



(4.5) 



and the corresponding Nemytskii operators F and B are obtained from ( |2.8| l and clearly satisfy Assump- 
tion pj] (if the domain of / is restricted to M+) and Assumption 2.4 (see |j40] Section 4]) respectively, 



where ( |4.2| i is used in the noise representation ( 1.2 1. The linear operator is given by 

A = V-DV(.)-V-q(.). 



(4.6) 



For a heterogeneous medium we used three parallel high permeability streaks. This could represent for 
example a highly idealized fracture pattern. We obtain the Darcy velocity field q by solving the system 



V-q = 
At 



(4.7) 



with Dirichlet boundary conditions = {0, 1 } x [0, 1] and Neumann boundary F^ = (0, 1 ) x {0, 1 } such 



that 



and 



m 
in 



{0}x[0,l] 
{l}x[0,l] 



F^ 



where p is the pressure, pL is dynamical viscosity and k the permeability of the porous medium. We have 
assumed that rock and fluids are incompressible and sources or sinks are absent, thus the equation 

r^(x) 



V q = V 



= 



comes from mass conservation. As in 
sentation ( |1.2| i 



, we take the following values for {qij}^ 



r>0. 



(4.8) 
in the repre- 

(4.9) 



Note that to have a trace class noise we need r > 2. In our simulation we use r = 2.01. To deal with 
high Peclet flows we discretize in space using finite volumes. We can write the semi-discrete finite volume 
discretization of (|4.31l as 



dX'' = {AhX'' +PhF{X^))+PhB{X^)dW, (4.10) 
(see Il33ll23l ). Figure 2(a) shows the convergence of SETDMO, SETDMl and semi-implicit schemes for 



the homogeneous porous medium. The scheme SETDMl seems to be more accurate for large time steps 
but for large time steps it has the same order of accuracy as the semi-implicit and SETDMO schemes. The 
temporal order is 0.54 for SETDMl scheme, 0.58 for SETDMO scheme and the semi-implicit scheme. We 
used 200 realizations and the convergence order is close to the 0.5, the predicted order of convergence in 
Theorem 2.7 A sample the 'true solution' is shown in Figure 2(b) with At — 1/1600 while the mean of the 



"true solution" for 200 realizations is shown in Figure [2(c)l 
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(a) 



(b) 



(c) 



Figure 2: (a) Convergence of the root mean square norm at T = 1 as a function of At with 200 realiza- 
tions with Ax = Ay = 1/150, Xo = for the homogeneous medium. The noise is white in time with ( |4.9| l 
and r = 2.01. The temporal orders of convergence in time are 0.54 for SETDMl and 0.58 for SETDMO 
semi-implicit schemes. In (b) we plot a sample of the 'true solution' for r — 2.01 with Af = 1 /1600 while 
(c) shows the mean of the true solution for 200 realizations. 



Figure [3(a)] shows the convergence of SETDMO and SETDMl schemes for the heterogeneous porous 
medium. It also shows that SETDMl is more accurate than SETDMO scheme for high time step size. The 



observed temporal order is 0.54 for SETDMl scheme and 0.58 for SETDMO scheme. Figure 3(c) shows 



the streamline of the velocity field. A sample the "true solution" is shown in Figure 3(b) with At ~ 1 /1600 



while the mean of the "true solution" for 200 realizations is shown in Figure 3(d) 



To conclude we have proved the errors estimates for the exponential based integrators and observed the 
predicted rate of convergence in the simulations. 
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(c) (d) 



Figure 3: (a) Convergence of the root mean square L? norm at T = 1 as a function of At with 200 real- 
izations with Ax = Ay = 1/150, Xq = for the heterogeneous medium. The noise is white in time with 
( |4.9| l and r = 2.01. The temporal orders of convergence in time are 0.54 for SETDMl scheme and 0.58 
for SETDMO. In (b) we plot a sample 'true solution' for r = 2.01 with At = 1/1600. In (c) we plot the 
streamline of the velocity field while (d) shows the mean of the "true solution" for 200 realizations. 
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